As noted, there are two sets of scaling factors that are used in these routines. The temporary scaling factors arguments (tmp_scale) are used exclusively to compensate for the rescaling for dimensional consistency. By contrast, although the net scaling factors (unscale or scale) had previously been used for the dimensional consistency factors in most cases, they are also used to change the units of variables - for example changing temperatures in to heat capacities or salinities into total salt contents. As a result, there are legitimate reasons to use both simultaneously.

Perhaps the documentation of the internal units needs to be cleaned up to describe all the various options, but I am still convinced that the code as revised is functionally correct.

As I am working through the comments on this PR, it is becoming increasingly clear that correctly rescaling the units is going to be very confusing if both the scale and unscale arguments are present in the various global_..._mean() and global_..._integral() routines. I am coming around to the suggestion that there should be an error issued if someone tries to use both simultaneously. I do not believe that we actually have any instances yet where both are used, and one can always do something like a conversion from integrated temperature to heat by multiplication after the fact or via a conversion argument to a register_diag_field() call.

I have now worked out how to very explicitly document all of the units of variables in the various spatial mean routines when both the unscale and tmp_scale arguments are present, and modified the comments in this commit accordingly, which will hopefully address many of your specific comments.

I have also corrected the double-counting bugs that you spotted, @herrwang0 .

I would appreciate it if you could take a look at the revised code and let me know if it addresses all of the points of confusion that you noted, or at least does as much as we can in a situation with so many combinations of options.

The two functions where the areal unscaling is conditionally included in scalefac when tmp_scale is absent are global_area_integral() and global_mass_integral(), in which case it is not just the variable that is being integrated that is being rescaled with the call, but also the area or mass. Although this pattern is inconsistent with the other functions, it is both functionally important and correct.

Thank you very much for your very thorough reviews, @herrwang0. I have now updated the comments in the code, several of them using your tables that were directly inspired by your comments. I am confident that the potentially confusing units of the variables in this PR are now as thoroughly documented as could be considered reasonable. I hope that you agree with me that all of your comments have been addressed and that this PR is now ready to go.

So this is perhaps not well documented in the comments, but without its unscale parameter, reproducing_sum() can only work properly with unscaled variables, so without the unscale or tmp_scale being passed to the routines in MOM_spatial_means.F90, the var argument must already be in arbitrary unscaled units (i.e., [a] and not [A ~> a]). With the last attempt at this PR, I had added [a] as possible units for var on line 163, but it was not explicit enough. In other words, the units in all the tables are correct, but the comments describing the var arguments are not explicit enough.

I have now revised this PR again to clarify this point.

So this is perhaps not well documented in the comments, but without its unscale parameter, reproducing_sum() can only work properly with unscaled variables, so without the unscale or tmp_scale being passed to the routines in MOM_spatial_means.F90, the var argument must already be in arbitrary unscaled units (i.e., [a] and not [A ~> a]). With the last attempt at this PR, I had added [a] as possible units for var on line 163, but it was not explicit enough. In other words, the units in all the tables are correct, but the comments describing the var arguments are not explicit enough.

I have now revised this PR again to clarify this point.

Thanks for the clarification! This PR looks good to me now.

And thanks for bearing with me for this unexpectedly-exploded review of an otherwise simple and straightforward PR. I appreciate the effort for making the documentation clearer, which I think also makes it easier to identify potential bugs.

Thank you very much for your thorough reviews, @herrwang0. The documentation of the code has clearly been significantly improved by this rigorous review process.